function  [logl,eu,nuhat] = tloglik4(para,data) 

global mprior psel 

if mprior == 4 
% parameter loading 

alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a = para(8);
rho_f = para(9);
rho_i  = para(10);
sd_a1  = para(11);
sd_f1  = para(12);
sd_i1   = para(13);
sd_a2  = para(14);
sd_f2  = para(15);
sd_i2   = para(16);
piess  = para(17);
gy     = para(18);
lnA_0  = para(19);
yst    = para(20);
p11 = para(21);
p22 = para(22);
p11_f = para(23);
p22_f = para(24);

bound  = zeros(length(para),2);
bound(8:10,1) = 1e-4; 
bound(19:20,1) = -100; 
bound(:,2) = [100;100;100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100;100;100;1;1;1;1]; 
     
crit1 = para-bound(:,1); 
crit2 = bound(:,2)-para; 

        if min(crit1)<0 || min(crit2)<0 
logl = -10000000; 
           return
        else


P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);

iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

eu =0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu =1; 
end

if eu ==0 
%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------
rho     = diag([rho_a;rho_f;rho_i]);

sd1       = [sd_a1;sd_f1;sd_i1];
sd2       = [sd_a2;sd_f2;sd_i2]; 
sd3 =sd1; sd4 =sd2; 


B1 = zeros(1,3); B2 = zeros(1,3);
B3 = zeros(1,3); B4 = zeros(1,3); 

B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i1a, i1f, i1i];
B3(1,:) =  [i2a, i2f, i2i];
B4(1,:) =  [i2a, i2f, i2i];

% Augment the state vector 
rho1 = diag([diag(rho);1]);
rho2 = diag([diag(rho);1]);
rho3 = diag([diag(rho);1]);
rho4 = diag([diag(rho);1]);
rho1(4,1) = rho_a; 
rho2(4,1) = rho_a; 
rho3(4,1) = rho_a; 
rho4(4,1) = rho_a; 
TC  = [0;0;0;gy]; 
sd1  = [sd1;sd_a1]; sd2  = [sd2;sd_a2];  
sd3  = [sd3;sd_a1]; sd4  = [sd4;sd_a2]; 
Vs1 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs2 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);
Vs3 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs4 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac2 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac3 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 
mac4 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 

% Markov Chain transition probability 
trans = P';  

% compute the stationary distribution 
A = [eye(length(trans))-trans;ones(1,length(trans))]; 
probs = inv(A'*A)*(A'*[zeros(length(trans),1);1]); 
Vs  = probs(1)*Vs1+probs(2)*Vs2+probs(3)*Vs3+probs(4)*Vs4;  

% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 2;              
nxmax = 4^nhist; 

Az      = [ yst ; piess ; iss]; 
Az1 = Az; Az2 = Az; Az3 = Az; Az4 = Az; 
Bz1     = [ mac1 ; B1(1,:) 0];    
Bz2     = [ mac2 ; B2(1,:) 0]; 
Bz3     = [ mac3 ; B3(1,:) 0];    
Bz4     = [ mac4 ; B4(1,:) 0]; 
Cz      = diag([0;0;0]); 

% condition on the first observation : assume the first regime 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz1*(TC+rho1*s00)); 
Ft        = Bz1*Vs*Bz1'; 
iFt         = Ft\eye(3);
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz1')*iFt*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz1'*iFt*(Vs*Bz1')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(4*nx,1);      
s0t      = zeros(4*nx,ns);    
ltrvecP0 = zeros(4*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat1= TC+rho1*s1t(i,:)';
   alphahat2= TC+rho2*s1t(i,:)';
   alphahat3= TC+rho3*s1t(i,:)';
   alphahat4= TC+rho4*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2; 
   Phat2    = rho2*P1t*rho2'+(diag(sd2))^2; 
   Phat3    = rho3*P1t*rho3'+(diag(sd3))^2; 
   Phat4    = rho4*P1t*rho4'+(diag(sd4))^2; 
   
   yhat1    = Az1 + Bz1*alphahat1;
   yhat2    = Az2 + Bz2*alphahat2; 
   yhat3    = Az3 + Bz3*alphahat3; 
   yhat4    = Az4 + Bz4*alphahat4; 
   
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   nu3      = data(t,:)-yhat3'; 
   nu4      = data(t,:)-yhat4'; 
   
   Ft1      = Bz1*Phat1*Bz1';
   Ft2      = Bz2*Phat2*Bz2'; 
   Ft3      = Bz3*Phat3*Bz3'; 
   Ft4      = Bz4*Phat4*Bz4'; 
   
   iFt1    = Ft1\eye(3);  iFt2=Ft2\eye(3); iFt3=Ft3\eye(3); iFt4=Ft4\eye(3); 
   
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2+probx2(2*nx+i)*yhat3+probx2(3*nx+i)*yhat4;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2')+probx2(2*nx+i)*(Ft3+yhat3*yhat3')+probx2(3*nx+1)*(Ft4+yhat4*yhat4'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*iFt1*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*iFt2*nu2'; 
   llik_ht(2*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft3))-0.5*nu3*iFt3*nu3'; 
   llik_ht(3*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft4))-0.5*nu4*iFt4*nu4'; 
    
   % updating step 
   s0t(i,:)    = (alphahat1+Phat1*Bz1'*iFt1*nu1')';
   s0t(nx+i,:) = (alphahat2+Phat2*Bz2'*iFt2*nu2')'; 
   s0t(2*nx+i,:) = (alphahat3+Phat3*Bz3'*iFt3*nu3')'; 
   s0t(3*nx+i,:) = (alphahat4+Phat4*Bz4'*iFt4*nu4')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*iFt1*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*iFt2*Bz2*Phat2'))'; 
   ltrvecP0(2*nx+i,:) = (ltrvec(Phat3-Phat3*Bz3'*iFt3*Bz3*Phat3'))'; 
   ltrvecP0(3*nx+i,:) = (ltrvec(Phat4-Phat4*Bz4'*iFt4*Bz4*Phat4'))'; 
    
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx));sum(probx2(2*nx+1:3*nx));sum(probx2(3*nx+1:4*nx))]; 
   
   % collapse the state vector : Kim and Nelson (1999)'s approximation 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (4*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat       = kron(eye(nx),[1 1 1 1]); 
   sumzeroprob  = selmat*zeroprob;
   del          = sumzeroprob==4;  
   nx4          = sum(del); 
   
   s1t = zeros(nx-nx4,ns);
   ltrvecP1 = zeros(nx-nx4,ntr); 
   probx = zeros(nx-nx4,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*4+1; 
       
       if sumzeroprob(i)<4 ; 
          % eliminate histories if all the subhistories have insignificant prob.    
          % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1)+probx2(i2+2)+probx2(i2+3); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:)+(probx2(i2+2)/probx(j))*s0t(i2+2,:)+(probx2(i2+3)/probx(j))*s0t(i2+3,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) +...
                          (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))+...
                          (probx2(i2+2))/(probx(j))*(ltrmat(ltrvecP0(i2+2,:)')+s0t(i2+2,:)'*s0t(i2+2,:))+... 
                          (probx2(i2+3))/(probx(j))*(ltrmat(ltrvecP0(i2+3,:)')+s0t(i2+3,:)'*s0t(i2+3,:))+...
                           -s1t(j,:)'*s1t(j,:))';
            j=j+1;
        
       end
       i=i+1;
   end
   end

   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
t=t+1;
end

logl = sum(llik_c);
%logl = llik_c; 
elseif eu>0  
logl = -1000000;
 end                  % determinacy case 
end                   %bound 

elseif mprior == 5 

% parameter loading 

alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a = para(8);
rho_f = para(9);
rho_i  = para(10);
sd_a1  = para(11);
sd_f1  = para(12);
sd_i1   = para(13);
sd_a2  = para(14);
sd_f2  = para(15);
sd_i2   = para(16);
gy     = para(17);
lnA_0  = para(18);
yst    = para(19);
p11 = para(20);
p22 = para(21);
p11_a = para(22);
p22_a = para(23);
lnP_0 = para(24);
sd_p1 = para(25);


bound  = zeros(length(para),2);
bound(18:19,1) = -100; 
bound(:,2) = [100;100;100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100;100;1;1;1;1;100;100]; 
     
crit1 = para-bound(:,1); 
crit2 = bound(:,2)-para; 

        if min(crit1)<0 || min(crit2)<0 
logl = -10000000; 
           return
        else


P_m  = [p11 1-p11;
       1-p22 p22];
   
P_a = [p11_a 1-p11_a;
       1-p22_a p22_a];   

P = kron(P_m ,P_a);



iss    = gy-log(beta1);

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

eu =0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu =1; 
end

if eu ==0 
%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];

 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------
rho1     = diag([rho_a;rho_f;rho_i;1;1]);
rho1(4,1)= rho_a; 


TC  = [0;0;0;gy;0]; 

sd1 = [sd_a1;sd_f1;sd_i1]; 
sd1  = [sd1;sd_a1;sd_p1]; sd2=[sd_a2;sd_f2;sd_i2;sd_a2;sd_p1]; 
sd3  = sd1;  sd4=sd2; 

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1 0; pi1a pi1f pi1i 0 1;i1a i1f i1i 0 1]; 
mac2 = [y2a y2f y2i 1 0; pi2a pi2f pi2i 0 1;i2a i2f i2i 0 1]; 

% Markov Chain transition probability 
trans = P';  

% compute the stationary distribution 
A = [eye(length(trans))-trans;ones(1,length(trans))]; 
probs = inv(A'*A)*(A'*[zeros(length(trans),1);1]); 
Vs1 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0;0]);
Vs2 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0;0]);
Vs3 = Vs1; Vs4=Vs2; 


Vs =probs(1)*Vs1+probs(2)*Vs2+probs(3)*Vs3+probs(4)*Vs4; 


% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 2;              
nxmax = 4^nhist; 

Az      = [ yst ; 0 ; iss]; 
Bz1     = mac1;     Bz2     = mac2; 

% condition on the first observation : assume the first regime w.l.o.g. 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 
s00(5)   = lnP_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu1      = data(1,:)'-(Az+Bz1*(TC+rho1*s00));
Ft       = Bz1*Vs*Bz1'; 

llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu1'*inv(Ft)*nu1;

s00      = s00+(Vs*Bz1')*inv(Ft)*nu1; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz1'*inv(Ft)*(Vs*Bz1')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(4*nx,1);      
s0t      = zeros(4*nx,ns);    
ltrvecP0 = zeros(4*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat = TC+rho1*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2;
   Phat2    = rho1*P1t*rho1'+(diag(sd2))^2;
   Phat3    = rho1*P1t*rho1'+(diag(sd3))^2; 
   Phat4    = rho1*P1t*rho1'+(diag(sd4))^2; 
    
   yhat1    = Az  + Bz1*alphahat;
   yhat2    = Az  + Bz2*alphahat; 
      
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
    
   Ft1      = Bz1*Phat1*Bz1'; Ft2 = Bz1*Phat2*Bz1'; Ft3 = Bz2*Phat3*Bz2';  Ft4 = Bz2*Phat4*Bz2';  
  
   yhat     = yhat + (probx2(i)+probx2(nx+i))*yhat1+(probx2(2*nx+i)+probx2(3*nx+i))*yhat2;
  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat1*yhat1')+probx2(2*nx+i)*(Ft3+yhat2*yhat2')+probx2(3*nx+1)*(Ft4+yhat2*yhat2');

   % compute likelihood 
   
   llik_ht(i,1)      = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1)   = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu1*inv(Ft2)*nu1'; 
   llik_ht(2*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft3))-0.5*nu2*inv(Ft3)*nu2'; 
   llik_ht(3*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft4))-0.5*nu2*inv(Ft4)*nu2'; 
  
 
   % updating step 
   s0t(i,:)    = (alphahat+Phat1*Bz1'*inv(Ft1)*nu1')'; s0t(nx+i,:) = (alphahat+Phat2*Bz1'*inv(Ft2)*nu1')'; 
   s0t(2*nx+i,:)= (alphahat+Phat3*Bz2'*inv(Ft3)*nu2');  s0t(3*nx+i,:) = (alphahat+Phat4*Bz2'*inv(Ft4)*nu2'); 
  
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))';     ltrvecP0(nx+i,:)   = (ltrvec(Phat2-Phat2*Bz1'*inv(Ft2)*Bz1*Phat2'))'; 
   ltrvecP0(2*nx+i,:)=(ltrvec(Phat3-Phat3*Bz2'*inv(Ft3)*Bz2*Phat3'))';  ltrvecP0(3*nx+i,:) = (ltrvec(Phat4-Phat4*Bz2'*inv(Ft4)*Bz2*Phat4'))';   
  
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx));sum(probx2(2*nx+1:3*nx));sum(probx2(3*nx+1:4*nx))];   
    % collapse the state vector : Kim and Nelson (1999)'s approximation 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (4*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat       = kron(eye(nx),[1 1 1 1]); 
   sumzeroprob = selmat*zeroprob;    
   del          = sumzeroprob==4;  
   nx4          = sum(del); 
   
   s1t = zeros(nx-nx4,ns);
   ltrvecP1 = zeros(nx-nx4,ntr); 
   probx = zeros(nx-nx4,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*4+1; 
       
       if sumzeroprob(i)<4 ; 
          % eliminate histories if all the subhistories have insignificant prob.    
          % collapse density 
            probx(j) = sum(probx2(i2:i2+3)); 
  s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:)+(probx2(i2+2)/probx(j))*s0t(i2+2,:)+(probx2(i2+3)/probx(j))*s0t(i2+3,:);
           ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) +...
                          (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))+...
                          (probx2(i2+2))/(probx(j))*(ltrmat(ltrvecP0(i2+2,:)')+s0t(i2+2,:)'*s0t(i2+2,:))+... 
                          (probx2(i2+3))/(probx(j))*(ltrmat(ltrvecP0(i2+3,:)')+s0t(i2+3,:)'*s0t(i2+3,:))-s1t(j,:)'*s1t(j,:))';
    
            j=j+1;
        
       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
t=t+1;
end
% logl=llik_c;
logl = sum(llik_c);
  if imag(logl)==0 
      logl=logl;
  else
      logl=-10000000;
  end 
elseif eu>0  
logl = -1000000;
 end                  % determinacy case 
end                   %bound 

end   % mprior 
   